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Abstract 

Two-phase  flow  dominated  by  capillary  effects  in  model  fibrous  media  is  studied  combining  pore-network  simulations  and  visualisations  on 
transparent  micromodels.  It  is  shown  that  the  process  of  liquid  water  invasion  in  a  hydrophobic  medium  can  be  simulated  using  the  classical  invasion 
percolation  algorithm  provided  that  the  contact  angle  (measured  in  air,  which  is  the  wetting  phase)  is  sufficiently  far  below  90°.  For  contact  angles 
approaching  90°,  changes  in  the  interface  local  growth  mechanisms  lead  to  changes  in  the  invasion  pattern. 

Then  it  is  shown  that  the  invasion  pattern  is  dramatically  different  in  a  hydrophilic  medium.  Impact  of  wettability  (hydrophobic  vs.  hydrophilic) 
on  evaporation  pattern  is  also  analysed. 

In  a  last  part,  implications  of  the  study  findings  on  the  water  management  problem  in  the  gas  diffusion  layers  (GDLs)  of  PEMFC  are  discussed. 
Our  results  provide  pore- scale  explanations  to  the  advantages  of  hydrophobic  GDLs. 

©  2007  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

PEM  fuel  cells  are  undergoing  intense  development.  Many 
groups  work  throughout  the  world  in  order  to  improve  their 
performances,  efficiency,  reliability,  manufacturability  and  cost- 
effectiveness.  In  this  context,  one  popular  way  for  studying  fuel 
cells  is  to  develop  computational  models  representing  the  com¬ 
plex  multi-physics  transport  processes  occurring  in  fuel  cells. 
This  leads  to  rather  complex  coupled  models  taking  into  account 
the  fluidic,  ionic,  electronic  and  thermal  transports  in  concert 
with  electrochemical  reactions,  e.g.  Refs.  [1-4]  and  references 
therein. 

While  much  progress  has  been  made  in  recent  years,  the 
development  of  truly  predictive  models  remains  a  challenge  due 
to  fundamental  modelling  issues.  Among  the  identified  issues, 
see  for  instance  the  discussion  in  Ref.  [5],  one  of  the  major  con¬ 
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cerns  is  the  modelling,  and,  more  generally,  the  understanding, 
of  two-phase  flow  in  the  gas  diffusion  layer  (GDL). 

The  GDLs  consist  of  an  anisotropic  fibrous  structure,  either 
in  the  form  of  a  random  structure  (carbon  paper)  or  woven  car¬ 
bon  fabrics  or  cloths.  The  thickness  of  a  GDL  typically  varies 
between  170  and  400  ptm  [6],  whereas  the  fibre  diameter  is 
typically  of  the  order  of  10  p,m  (carbon-paper  GDL)  leading 
to  a  distribution  of  pore  sizes  ranging  from  a  few  microns  to 
tens  of  microns  [7].  The  overall  porosity  of  a  GDL  is  in  the 
range  70-80%.  Thus  a  GDL  appears  as  a  highly  anisotropic 
fibrous  structure  very  thin  in  the  main  direction  of  transport  with 
thickness  of  only  a  few  tens  of  fibre  diameters.  Depending  on 
the  applications  considered,  one  or  several  microporous  layers 
(MPLs)  may  be  attached  to  the  fibrous  structure  [6],  in  order  to 
improve  the  electrical  contact  with  the  catalyst  layer  (CL)  and 
the  control  of  the  water  management.  However,  the  microstruc¬ 
ture  of  the  MPL  is  quite  different  from  that  of  the  macro-porous 
substrate  described  above  and  the  specific  effect  of  the  MPL  is 
not  studied  in  this  paper. 

The  GDL  are  commonly  teflonized,  rendering  them 
hydrophobic.  The  degree  of  hydrophobicity  is  a  priori  depen¬ 
dent  upon  the  amount  of  Teflon  added  to  the  GDL.  A  classic 
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Fig.  1.  Hydrophilic  (i 9<jt/2 )  and  hydrophobic  ( 6>tt/2 )  surfaces. 


way  of  characterizing  the  hydrophobicity  is  to  use  the  concept 
of  contact  angle,  which  is  the  angle  6  between  the  solid  phase 
and  the  interface  between  the  two  fluids  (see  Fig.  1).  When 
the  two  fluids  are  liquid  water  and  air,  the  surface  is  said  to 
be  hydrophilic  when  this  angle  measured  in  the  water  is  lower 
than  90°  and  hydrophobic  when  0>9O°.  The  contact  angle  for 
water  on  a  Teflon  flat  surface  is  roughly  110°,  whereas  water 
on  a  carbon  surface  is  close  to  80°  (thus  hydrophilic).  The  mea¬ 
surement  of  contact  angle  using  the  sessile  drop  method  on  a 
teflonized  GDL  carbon  cloth  is  reported  in  Ref.  [8].  It  is  found 
to  be  133°.  The  increased  contact  angle  as  compared  to  a  flat 
surface  is  attributed  to  the  surface  roughness  of  the  GDL.  The 
contact  angle  is  generally  measured  at  room  temperature.  Vari¬ 
ations  of  contact  angle  with  temperature  are  a  priori  expected. 
It  would  be  therefore  interesting  to  measure  it  at  temperatures 
representative  of  PEM  fuel  cell  operating  conditions  (~80°C). 
While  a  major  change  in  wettability  (i.e.,  from  hydrophobicity  to 
hydrophilicity)  is  not  expected,  it  would  be  interesting  to  evalu¬ 
ate  the  influence  of  varying  contact  angles.  Another  issue  is  that 
the  contact  angle  exhibits  hysteresis,  meaning  that  the  advanc¬ 
ing  (0a)  and  receding  (0r)  contact  angles  may  be  very  different, 
with  the  advancing  contact  angle  being  typically  greater  than  the 
receding  one.  The  difference  is  on  the  order  of  14°  for  distilled 
water  on  a  Teflon  flat  surface  according  to  the  measurements 
reported  in  Ref.  [8]  (with  0a  ~  127°  and  0r  ~  113°).  As  pointed 
out  in  Ref.  [9],  measurements  of  contact  angle  on  a  pure  Teflon 
flat  plate  or  at  the  external  surface  of  a  GDL  are,  however,  not 
necessarily  representative  of  the  internal  contact  angle,  i.e.  the 
average  contact  angle  within  the  pore  space  of  a  GDL.  The  val¬ 
ues  of  the  internal  contact  angle  in  GDL  materials  reported  in 
Ref.  [9]  are  in  the  range  [88°,  101°],  depending  on  the  GDL 
materials  considered  and  therefore  significantly  lower  than  the 
expected  value  on  a  pure  Teflon  plate.  As  we  shall  see,  the  wet¬ 
tability  (i.e.,  the  hydrophobic  or  hydrophilic  nature)  of  the  GDL 
has  a  great  impact  on  the  water  transport,  especially  in  the  range 
of  contact  angles  close  to  90°.  The  motivation  for  rendering  the 
GDL  hydrophobic  is  generally  associated  with  the  idea  that  the 
hydrophobicity  will  limit  the  flooding  of  the  GDL,  for  exam¬ 
ple  by  forcing  the  water  to  agglomerate  at  the  free  surface  of 
the  GDL  [5].  The  results  presented  in  this  paper  will  shed  some 
light  on  these  different  aspects. 

Another  important  issue  is  the  general  validity  of  macro¬ 
scale  two-phase  flow  models  commonly  used  when  modelling 
transport  phenomena  in  GDLs  and/or  between  the  GDL  and 
the  catalyst  layer,  e.g.  Ref.  [10].  Generally,  these  models  are 
the  conventional  models  developed  for  studying  flow  in  “usual” 
porous  media,  such  as  soils  or  sedimentary  rocks.  They  are 
based  on  the  so-called  generalized  Darcy’s  law  and  the  concepts 
of  macroscopic  capillary  pressure  and  relative  permeability, 


see  for  instance  Ref.  [11].  In  principle,  this  type  of  model  is 
intimately  related  to  the  concept  of  representative  elementary 
volume  (REV),  e.g.  Ref.  [12].  This  volume  is  supposed  to  be 
large  enough  for  being  representative  of  the  micro  structure  but 
also  sufficiently  small  compared  to  the  porous  domain  size.  As 
mentioned  before,  the  thickness  of  a  GDL  is  typically  on  the 
order  of  a  few  tens  of  fibre  or  pore  sizes,  which  is  presumably 
on  the  order  or  less  than  the  REV  size.  This  raises  fundamental 
questions  regarding  the  meaning  to  be  given  to  the  fields,  such  as 
saturation,  pressures  and  velocities,  computed  across  the  GDL 
using  this  type  of  model.  Hence,  using  this  type  of  modelling 
is  not  only  difficult  because  of  the  lack  of  data  for  the  various 
fibrous  materials  forming  a  GDL,  as  pointed  out  for  example 
in  Ref.  [5],  but  raises  also  a  more  fundamental  question  regard¬ 
ing  the  validity  of  the  continuum  approach  applied  to  ultra-thin 
porous  media  such  as  the  GDL. 

The  fact  that  this  approach  has  not  yet  been  validated  at  all  for 
the  GDLs  is  a  well-identified  problem,  e.g.  Refs.  [4,5,7].  Hence, 
there  is  a  need  for  a  better  understanding  of  two-phase  flows 
in  GDLs,  in  particular  from  a  micro-scale  perspective.  In  this 
respect,  visualisations  of  liquid  water  in  a  GDL  and  of  its  evolu¬ 
tion  are  certainly  useful.  Although  in  situ  visualisations,  e.g.  Ref. 
[13],  are  certainly  the  best  option  to  truly  understand  the  physics 
of  two-phase  flows  in  GDLs  of  PEMFC,  ex  situ  visualisations, 
e.g.  Ref.  [14],  are  also  useful  for  improving  our  understanding  of 
two-phase  flows  and  providing  data  for  comparison  with  mod¬ 
els.  On  the  modelling  side,  one  option  is  to  develop  simulations 
directly  at  the  scale  of  the  microstructure,  see  for  example  Ref. 
[15].  In  addition  to  a  better  understanding  of  the  physics,  this 
type  of  simulation  can  be  useful  for  evaluating  the  validity  of 
the  traditional  continuum  approach  and  also  for  computing  data 
difficult  to  measure  in  very  thin  systems,  such  as  the  capillary 
pressure-saturation  curve  or  the  relative  permeabilities.  In  this 
context,  although  limited  to  2D  model  porous  media,  the  results 
presented  in  this  paper  can  be  seen  as  a  contribution  to  a  bet¬ 
ter  understanding  of  two-phase  flow  in  model  fibrous  media  in 
relation  with  the  general  problems  of  two-phase  flow  modelling 
in  GDLs  and  the  impact  of  wettability  mentioned  above. 

The  paper  is  organized  as  follows:  first  we  present  in  Section 
2  a  model  of  liquid  water  quasi- static  invasion  in  a  hydrophobic 
model  fibrous  medium.  Results  of  simulation  are  compared  to 
experimental  visualisations.  In  Section  3,  we  concentrate  on  the 
influence  of  contact  angle  on  the  invasion  pattern.  In  Section  4, 
we  recall  recent  results  regarding  the  impact  of  wettability  on 
evaporation  patterns.  Impacts  of  the  results  on  two-phase  flows 
in  GDLs  are  discussed  in  Section  5. 

2.  Quasi-static  invasion  of  liquid  water  in  a  model  2D 
fibrous  medium 

As  a  matter  of  fact,  how  water  flooding  occurs  in  a  GDL  is 
not  perfectly  clear.  For  example,  Nam  and  Kaviany  [16],  tend 
to  consider  that  liquid  water  invades  the  GDL  as  a  result  of  a 
condensation  process.  In  other  studies,  e.g.  Ref.  [17],  the  conden¬ 
sation  process  is  not  taken  into  account  and  water  enters  directly 
as  a  liquid  phase  into  the  GDL  from  the  CL/GDL  interface.  The 
invasion  in  liquid  phase  from  one  face  of  the  fibrous  material 
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is  also  the  situation  considered  in  the  pore-scale  simulations  of 
Schulz  et  al.  [15],  and  the  visualisation  experiment  presented  in 
Ref.  [14].  Although  the  condensation  process  and  liquid  phase 
invasion  from  the  CL  are  likely  to  occur  simultaneously  in  an 
operating  fuel  cell,  we  consider  in  this  paper  the  same  invasion 
scenario  as  in  Refs.  [14]  or  [17]  and  therefore  suppose  that  liquid 
water  enters  the  fibrous  domain  from  one  side.  This  situation  is 
sketched  in  Fig.  2.  Whatever  the  exact  scenario  in  an  operating 
fuel  cell,  this  situation  is  interesting,  in  any  case,  for  evaluat¬ 
ing  two-phase  flow  models.  Initially,  the  fibrous  material  is  dry 
with  air  occupying  the  pore  space.  Then  liquid  water  is  injected 
into  the  pore  space  from  one  side  of  the  porous  domain  until 
breakthrough,  i.e.  until  water  reaches  the  side  opposite  to  the 
injection  side  (which  corresponds  to  the  bipolar  plate  channel  in 
a  PEMFC). 

In  this  section,  we  assume  that  the  GDL  is  hydrophobic. 
Hence  liquid  water  is  the  non- wetting  fluid  and  air  the  wetting 
one.  The  process  of  immiscible  displacement  of  a  wetting  fluid 
by  a  non-wetting  one  is  usually  called  drainage,  e.g.  Ref.  [11]. 
Therefore  the  invasion  of  liquid  water  into  a  dry  hydrophobic 
fibrous  material  corresponds  to  a  drainage  process.  This  pro¬ 
cess  was  studied  in  detail  in  Ref.  [18].  As  shown  in  Ref.  [18], 
this  process  is  controlled  by  the  competition  between  capillary 
effects  and  viscous  effects  (assuming  negligible  gravity  effects). 
The  invasion  pattern  depends  on  two  parameters:  the  capillary 
number  Ca  and  the  ratio  of  dynamic  viscosities  of  the  two  flu¬ 
ids  M  =  /rnw//xw  ~  10-3  for  the  water/air  system.  The  capillary 
number  Ca ,  which  is  the  ratio  between  viscous  forces  acting  at 
the  pore  scale  in  water  and  capillary  forces,  is  expressed  as 


Ay  cos  0W 

where  y  is  the  surface  tension,  /xnw  is  the  liquid  water  dynamic 
viscosity,  0W  is  the  contact  angle  taken  in  the  wetting  fluid  (air), 
A  is  the  cross-sectional  area  of  the  sample  (area  of  injected  face) 
and  Q\  is  the  injected  flow  rate. 


Fig.  3.  Model  fibrous  medium,  Voronoi  diagram  and  pore  network. 

To  estimate  Ca  in  the  case  of  PEMFC s,  we  assume  that  all 
the  water  produced  within  the  CL  invades  the  GDL.  This  leads 
to  express  Ca  as 

Ca  =  Mnw/MH2°  (2) 

ycos  0w2Fp\ 

where  F  is  the  Faraday’s  constant,  Mh2o  is  the  water  molec¬ 
ular  weight,  p\  is  the  liquid  water  density  and  i  is  the  current 
density.  With  /=lAcm-2  and  0W  =  70°,  Eq.  (2)  leads  to 
Ca  ~  10-7-10-8.  According  to  Ref.  [18],  the  displacement  for 
such  low  values  of  the  capillary  number  is  dominated  by  cap¬ 
illary  effects  and  leads  to  a  regime  called  capillary  fingering. 
In  this  low  capillary  number  limit,  the  drainage  process  can  be 
simulated  using  the  invasion  percolation  (IP)  algorithm  [19],  at 
least  when  the  contact  angle  0W  is  sufficiently  low  (see  Section 

3). 

In  the  rest  of  this  section  we  show  how  this  algorithm  can 
be  implemented  taking  as  an  example  the  invasion  within  a  2D 
model  fibrous  medium. 

2.7.  Model  fibrous  medium  and  associated  pore-network 
model 

As  depicted  in  Fig.  3,  the  system  under  consideration  is 
formed  by  an  array  of  equal- sized  non-overlapping  cylinders 
(disks)  randomly  distributed.  This  is  similar  to  a  cross-section 
through  an  anisotropic  arrangement  of  equal  diameter  fibres.  The 
disk  configuration  is  generated  using  the  Metropolis  algorithm 
as  described  in  Ref.  [20]. 

The  next  step  is  to  map  the  pore  space  onto  a  discretized  inter¬ 
connected  network  of  sites  (pores)  connected  by  bonds  (throats). 
To  this  end,  we  begin  with  the  computation  of  the  Voronoi  dia¬ 
gram,  following  [20].  This  leads  to  a  partition  of  the  domain  into 
regions  limited  by  polygons,  as  shown  in  Fig.  3.  Each  polygon 
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surrounds  a  disk  so  that  each  point  inside  the  polygon  is  closer 
to  this  disk  than  to  any  other  disk  in  the  system.  A  segment  of  a 
Voronoi  polygon  is  formed  by  points  that  are  at  the  same  distance 
from  two  disks.  The  polygon  segments  and  polygon  segments 
intersections  are  identified  as  the  bonds  and  sites  (pores)  of  the 
pore  network,  respectively.  The  size  of  a  bond  is  identified  as 
the  diameter  of  the  smallest  circle  with  center  on  the  bond  and 
in  contact  with  the  two  disks  adjoint  to  the  polygon  segment 
forming  the  bond.  In  other  terms,  if  Ci  and  C2  are  the  centers  of 
the  two  disks  adjoining  the  bond,  the  bond  size  i  is  expressed 
as  t  =  IC1C2I  —  D,  where  IC1C2I  and  D  are  the  Euclidian  dis¬ 
tance  and  disk  diameter,  respectively.  For  simplicity,  no  volume 
is  associated  with  bonds,  which  act  only  as  capillary  barriers,  as 
will  be  seen  below.  Volumes  are  associated  with  pores  and  corre¬ 
spond  to  the  regions  of  pore  space  limited  by  the  bond  minimum 
apertures. 


2.2.  Drainage  algorithm 


As  in  any  bond  invasion  percolation  model  on  a  network, 
the  first  step  is  to  assign  an  invasion  capillary  threshold  to  each 
bond.  In  our  case,  the  threshold  is  associated  with  the  aperture 
minimum  between  the  two  disks  limiting  the  bond.  According 
to  Laplace’s  law,  the  threshold  capillary  pressure  reads 


P\[ 1  —  2  "Y  cos  Qyy 


(3) 


where  e  is  the  length  of  the  cylinders.  Note  that  the  simulations 
are  performed  in  2D  (which  corresponds  to  e  =  00)  whereas  the 
experiment  (see  below)  is  carried  out  with  finite  equal  length 
cylinders  with  e>i.  Since  e>i  and  e  is  uniform,  the  capillary 
pressure  threshold  is  essentially  controlled  by  the  minimum  gap 
i  between  the  cylinders  also  in  the  experiments.  Hence,  we  can 
use  the  approximation  Pth  ~  (2 y  cos  0w)/l  both  in  the  simulation 
and  to  interpret  the  experiment. 

Implementing  the  IP  algorithm  is  now  straightforward.  Ini¬ 
tially,  all  bonds  are  occupied  by  the  wetting  fluid  (air).  The 
non- wetting  fluid  invades  the  network  from  one  of  its  four  sides; 
the  bottom  side  in  our  case  as  depicted  in  Fig.  3.  The  algorithm 
then  consists  of  repeating  two  steps: 


1.  Identify  the  bond  adjacent  to  the  invaded  region  that  has  the 
lowest  invasion  threshold. 

2.  Invade  the  identified  bond  and  fill  the  adjacent  pore  up  to  the 
next  cylinder  minimum  gaps. 


With  the  two  additional  rules: 


(3)  Bonds  not  belonging  to  the  percolation  cluster  skeleton  can¬ 
not  be  invaded. 

(4)  Bonds  belonging  to  trapped  wetting  fluid  bond  cluster  can¬ 
not  be  invaded. 


The  percolation  cluster  refers  to  the  set  of  bonds  occupied 
by  the  defender  (wetting)  phase  and  connected  to  the  exit  side 
through  at  least  one  path  of  defender  occupied  bonds  and  pores. 


The  invasion  stops  when  the  non-wetting  fluid  forms  a 
spanning  cluster  from  bottom  to  topsides.  No  flow  boundary 
conditions  are  imposed  on  lateral  sides  of  network  to  fit  in  with 
the  experimental  conditions.  The  percolation  cluster  is  identified 
using  a  depth  search  algorithm,  e.g.  Ref.  [21]. 

Results  of  drainage  simulation  based  on  this  algorithm  are 
compared  to  experimental  visualisations  in  a  machined  micro¬ 
model  after  the  next  section,  which  describes  the  experiments. 

2.3.  Experiments 

A  numerically  generated  hydrophobic  model  fibrous  medium 
has  been  constructed  using  a  digitally  controlled  milling 
machine.  The  fabrication  procedure  can  be  summarized  as  fol¬ 
lows,  see  Ref.  [22]  for  more  details.  First,  1  mm  diameter,  1  mm 
long  cylindrical  holes  are  drilled  into  a  Plexiglas  plate  with  a 
micromiller.  The  machining  precision  is  equal  to  10  p,m  in  every 
(v;  y;  z)  direction.  The  position  of  cylinders  is  generated  numeri¬ 
cally  using  the  Metropolis  algorithm  as  for  the  numerical  model. 
More  precisely,  the  same  distribution  of  cylinders  is  used  both 
in  the  experiment  and  the  simulation.  Note  that  the  micromodel 
1-mm  depth  is  lower  than  the  capillary  length  (which  is  on  the 
order  of  1.5  mm)  in  order  to  avoid  gravity  effects  on  the  liquid 
distribution  within  the  depth.  Thus,  apart  from  the  uncertainties 
introduced  in  the  course  of  fabrication,  the  pore  space  geometry 
in  the  model  plane  is  the  same  in  the  experiment  and  the  simu¬ 
lation.  Then  the  plate  is  machined  so  as  to  contain  the  negatives 
of  fluid  supply  and  exit  channels  along  the  cylinders  domain. 
The  next  step  is  to  mould  the  machined  Plexiglas  plate  with 
Rhodorsil®  RTV-2  (a  silicone  elastomer).  This  gives  a  RTV  slice 
containing  the  array  of  solid  cylinders  with  the  supply  and  exit 
channels.  Then  the  RTV  slice  is  sandwiched  between  a  metallic 
support  plate  and  a  Plexiglas  transparent  cover  plate.  Tightness 
is  obtained  by  compressing  a  little  the  RTV  between  the  metallic 
and  Plexiglas  plates  thanks  to  a  series  of  equally  spaced  and  con¬ 
trolled  clamping  screws.  The  exit  channel  along  the  topside  of 
the  cylinder  array  is  open  to  ambient  pressure,  while  the  entrance 
channel  along  the  bottom  side  is  connected  to  a  screw-driven 
pump.  The  micromodel  is  initially  fully  saturated  with  air  and 
put  on  a  horizontal  table  so  as  to  avoid  the  influence  of  gravity 
effects. 

Experiments  are  carried  out  by  slowly  injecting  water  into 
the  micromodel  from  the  bottom  side.  A  constant  temperature 
of  22  °C  is  maintained  during  the  experiment.  The  experiment  is 
recorded  with  a  CCD  camera  set  above  the  micromodel.  There 
are  some  unavoidable  experimental  uncertainties  concerning  the 
exact  reproduction  of  the  numerical  array  of  cylinder  coming 
from  the  machining  precision,  elastic  deformations  of  the  mate¬ 
rials,  the  plate’s  parallelism,  etc.  However,  as  we  shall  see,  this 
does  not  prevent  a  favourable  comparison  between  the  experi¬ 
mental  and  numerical  results. 

The  results  reported  in  this  paper  have  been  obtained  for  a 
2.5  cm  x  2.5  cm  model  containing  about  23  x  23  cylinders,  as 
shown  in  Fig.  4.  Note  that  the  contact  angle,  taken  in  the  water 
phase,  is  ~110°  for  RTV/water/air  [23],  which  is  comparable 
to  the  value  over  teflonized  surfaces.  Thus,  0W  ~  70°.  The  water 
injection  flow  rate  imposed  in  the  experiment  is  0.75  ml  h-1, 
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Fig.  4.  Top  view  of  physical  micromodel. 


which  leads  to  Ca~  5  x  10-7.  It  can  be  noticed  that  the  pore 
sizes  in  the  experiments  are  larger  than  in  a  GDL  material.  How¬ 
ever,  this  has  no  effect  on  the  pattern  provided  that  the  capillary 
number  is  sufficiently  small. 

2.4.  Invasion  patterns 

As  mentioned  before,  the  drainage  simulation  has  been  per¬ 
formed  using  the  IP  algorithm  described  in  Section  2.2  over  an 
array  of  disks  in  principle  identical  to  the  experimental  model 
described  in  Section  2.3. 

Fig.  5a  shows  the  experimental  invasion  patterns  obtained  at 
different  times  in  the  physical  model.  Fig.  5b  shows  the  numeri¬ 
cal  prediction  using  the  IP  algorithm  on  the  network  constructed 
from  the  Voronoi  diagram  at  comparable  saturations.  Overall,  the 
comparison  in  Fig.  3  shows  a  good  agreement  between  the  exper¬ 
imental  observations  and  the  numerical  computations.  Although 
the  comparison  could  seem  to  be  qualitative  only,  it  should  be 
pointed  out  that  direct  comparison  between  experimental  and 
numerical  results  for  the  same  microstructure  are  scarce  in  the 
literature.  Often,  the  comparison  is  performed  only  between  sys¬ 
tems  having  the  same  statistical  properties  (porosity,  pore  and 
bond  size  distributions,  average  coordinance,  etc.).  Thus,  it  is 
interesting  to  see  that  the  IP  algorithm  leads  to  good  results  on  a 
particular  realisation,  which  is  as  much  as  possible  identical  to 
the  experiment. 

Some  differences  can  be  observed,  however,  for  example  at 
the  beginning  of  the  drainage  process  or  at  breakthrough  (Fig.  5). 
These  differences  could  be  attributed  to  the  previously  men¬ 
tioned  uncertainties  induced  by  the  micromodel  fabrication  or 
to  viscous  effects.  The  viscous  effects  are  neglected  in  the  sim¬ 
ulations.  Although  the  experiment  is  conducted  at  a  very  low 
capillary  number,  a  slight  influence  of  viscous  effects  could 
affect  the  invasion  patterns  compared  to  the  purely  quasi- static 


ones  obtained  in  the  simulation,  see  Ref.  [18].  A  possible  source 
of  discrepancy  may  be  found  in  the  relatively  high  contact  angle 
characterizing  the  experimental  system.  This  issue  is  explored 
in  detail  in  the  next  section. 

Before,  we  note  that  the  capillary  fingering  regime  illustrated 
here  is  qualitatively  consistent  with  the  liquid  clusters  that  have 
been  observed  in  a  GDL  by  means  of  synchrotron  X-ray  radio- 
graphies  [13].  The  invasion  patterns  reported  in  Ref.  [14]  also 
resemble  IP  patterns.  However,  as  pointed  out  in  Ref.  [15],  the 
injected  flow  rate  considered  in  Ref.  [14],  is  about  two  orders 
of  magnitude  greater  than  the  flow  rate  expected  in  an  operating 
fluid  cells,  which  leads  to  a  capillary  number  of  the  order  of 
10-6-10-5.  Hence  the  invasion  process  considered  in  Ref.  [14] 
is  likely  to  be  affected  by  viscous  effects. 

3.  Impact  of  wettability  on  the  invasion  pattern 

An  important  issue  in  relation  with  the  crucial  problem  of 
water  management  in  PEMFCs  is  the  influence  of  wettabil¬ 
ity  conditions  on  invasion  pattern.  Since  the  capillary  number 
characterizing  the  water  invasion  in  GDLs  is  expected  to  be 
very  small,  we  consider  capillarity-dominated  regimes  only  (no 
viscous  effects).  In  this  quasi- static  limit,  we  investigate  the 
influence  of  contact  angle  on  the  invasion  pattern,  considering 
again  a  2D  model  fibrous  medium. 

To  facilitate  the  numerics,  this  model  is  slightly  different  from 
the  one  considered  in  the  previous  section.  As  depicted  in  Fig.  6, 
disks  (the  solid  phase  of  the  fibrous  medium)  are  placed  on  a 
square  lattice  of  constant  lattice  spacing  a.  Disks  radii  are  ran¬ 
domly  distributed  according  to  a  uniform  distribution  law  in  the 
range  [Rmm,  Rmaxh  For  the  simulations  presented  in  the  follow¬ 
ing,  we  considered  an  array  of  14  x  20  disks  with  Rmin  =  Rm ax/3 
and  Rmax  =  0.48a.  This  model  is  similar  to  the  ones  considered 
in  Ref.  [24]  but  adapted  for  flow  rate  controlled  invasion  rather 
than  pressure-controlled  invasion. 

The  interface  between  the  invading  and  displaced  fluids  con¬ 
sists  of  a  sequence  of  arcs  between  pairs  of  disks  (in  Fig.  6, 
which  shows  the  invading  fluid  at  the  entrance  of  the  system 
before  the  first  invasion,  the  invading  fluid  is  in  grey,  the  dis¬ 
placed  fluid  is  in  white  and  the  interface  is  visible  as  a  series  of 
arcs  joining  the  bottom  row  of  disks).  At  each  step  of  the  inva¬ 
sion,  one  “interfacial  pore”  of  the  system  is  invaded.  A  “pore” 
is  defined  as  the  region  of  the  pore  space  located  inside  the 
square  formed  by  the  segments  joining  two  neighbour  disks  on 
the  square  lattice.  An  interfacial  pore  is  a  pore  not  yet  invaded 
containing  at  least  one  arc.  Interfacial  pores  are  indicated  with 
a  “A”  in  Fig.  6.  To  determine  which  interfacial  pore  should  be 
invaded  at  a  given  step  of  the  invasion,  we  determine  the  “criti¬ 
cal”  invasion  curvature  radius  of  each  interfacial  arc.  As  sketched 
in  Fig.  7,  this  consists  in  studying  the  arc  curvature  radius  as  the 
arc  moves  between  the  two  disks  it  intercepts.  As  discussed  in 
Ref.  [24],  there  are  three  main  local  growth  mechanisms  con¬ 
trolling  the  invasion.  The  first,  as  illustrated  in  Fig.  7a,  is  called 
“burst”.  It  corresponds  to  the  burst  of  the  meniscus  when  the 
local  capillary  pressure  exceeds  the  maximum  local  capillary 
pressure  compatible  with  a  stable  arc  joining  the  pair  of  disks 
under  consideration,  see  Ref.  [24]  or  [22]  for  more  details.  We 
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(a)  Experiment 


(b)Simulation 


Fig.  5.  (a)  Displacement  patterns  observed  during  the  slow  drainage  experiment  within  the  model  porous  medium.  Water  (in  white)  injected  from  the  bottom  edge 
displaces  air  (in  dark  grey),  (b)  Simulated  displacement  patterns  using  the  IP  model  on  the  Voronoi  network.  Invading  phase  in  black. 
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Fig.  6.  Model  fibrous  medium  formed  by  a  regular  array  of  disks  of  variable 
diameter  (invading  fluid  in  grey,  displaced  fluid  in  white).  The  interface  is  visible 
as  a  series  of  arcs  joining  the  bottom  row  of  disks. 


note  rb  the  arc  curvature  radius  corresponding  to  the  last  arc 
stable  before  burst.  As  sketched  in  Fig.  7b-d,  the  local  capil¬ 
lary  pressure  corresponding  to  the  arc  burst  may  be  not  reached 
because  another  event  occurs  before:  the  considered  arc  touches 


a  third  disk  during  its  growth,  Fig.  7b,  or  coalesces  with  another 
arc,  Fig.  7c  and  d.  If  this  happens,  we  compute  the  arc  radius 
of  curvature  corresponding  to  each  of  these  local  events  with 
rt  the  radius  corresponding  to  the  touch  event  and  rc  the  radius 
corresponding  to  the  coalescence  event,  if  any.  These  radii  are 
computed  as  functions  of  contact  angle,  disks  radii  and  distances 
between  disks  [24,22] .  We  affect  to  each  interfacial  arc  the  radius 
ra  =  max(rb,  rt,  rc).  This  gives  a  hierarchy  of  radii.  At  each  step 
of  the  invasion,  the  interfacial  pore  associated  with  the  arc(s) 
corresponding  to  max  (ra),  i.e.  the  minimum  invasion  threshold 
local  capillary  pressure,  is  invaded.  The  corresponding  “unsta¬ 
ble”  arc(s)  is  (are)  suppressed  from  the  list  of  interfacial  arcs  and 
new  arcs  are  created  at  the  periphery  of  the  invaded  pore.  This 
procedure  is  repeated  until  breakthrough.  Again  the  interested 
reader  can  refer  to  [22]  or  [24]  for  more  details  on  the  invasion 
procedure. 

Fig.  8  shows  the  evolution  of  the  invasion  pattern  as  a  func¬ 
tion  of  the  contact  angle  Oe  (contact  angle  taken  in  the  displaced 
fluid).  If  liquid  water  is  supposed  to  be  the  invading  fluid, 
a  hydrophobic  fibrous  medium  corresponds  to  0<$e<90°, 
whereas  a  hydrophilic  one  corresponds  to  90°  <0e  <  180°.  As 
can  be  seen  from  Fig.  8,  the  wettability  has  a  major  impact  on  the 
invasion  pattern.  For  the  low  contact  angles  (<^very  hydropho- 


Fig.  7.  Growth  of  interfacial  arcs  during  the  invasion.  The  arrows  indicate  the  meniscus  displacement  direction.  Local  invasion  events:  (a)  burst,  (b)  touch,  and  (c) 
and  (d)  arc  coalescence. 


Fig.  8.  (A)  Evolution  of  invasion  pattern  at  breakthrough  in  a  model  fibrous  medium  as  a  function  of  contact  angle  in  the  displaced  fluid  6>e-  Invading  phase  in  grey, 
displaced  phase  in  white.  (B)  Detailed  evolution  of  the  pattern  in  the  range  of  contact  angles  [82°,  98°]. 
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Fig.  9.  Probability  (in  %)  of  interface  local  growth  mechanisms  as  a  function  of  contact  angle  in  the  displaced  fluid. 


bic  GDL),  capillary  fingering  is  obtained  as  in  Section  2.  For 
contact  angles  greater  than  about  100°  (<^hydrophilic  GDL),  a 
compact  pattern  characterized  by  an  almost  flat  invasion  front  is 
obtained.  It  is  interesting  to  mention  that  the  transition  between 
the  two  patterns  is  not  abrupt  at  90°  but  mainly  occurs  in  the 
range  of  contact  angles  [80°,  100°]. 

The  change  in  the  invasion  pattern  can  be  explained  by  the 
change  in  the  interface  local  growth  mechanisms.  We  have 
plotted  in  Fig.  9  the  probability  of  each  elementary  local 
growth  mechanisms  as  a  function  of  contact  angle  (for  the  sin¬ 
gle  realisation  of  the  medium  considered  here,  the  probability 
of  a  given  elementary  event  is  simply  njnv  where  nQ  is  the 
number  of  pores  invaded  as  a  result  of  the  considered  local 
event/mechanism  and  nv  is  the  number  of  pore  invasion  to  reach 
breakthrough). 

As  can  be  seen  from  Fig.  9,  burst  is  the  dominant  mechanism 
for  low  contact  angle,  and  the  invasion  can  be  simulated  using 
the  IP  algorithm  as  in  Section  2.  For  contact  angles  greater  than 
about  100°,  coalescence  of  menisci  is  the  dominant  mechanism 
and  the  invasion  can  be  simulated  very  simply  by  a  flat  travel¬ 
ling  front  (note  however  that  for  sufficiently  high  contact  angles 
above  100°,  film  effects  can  complicate  the  invasion  with  certain 
pore  geometries,  e.g.  Ref.  [25]  and  references  therein). 

On  the  basis  of  the  results  shown  in  Figs.  8  and  9,  the  slight 
differences  observed  between  the  experiment  and  the  simulation 
in  Section  2  can  be  attributed  to  the  fact  that  the  local  displace¬ 
ment  in  the  experiment  is  not  completely  dominated  by  the  burst 
mechanism  (as  opposed  to  the  simulation  in  Section  2  where 


only  the  burst  mechanism  is  taken  into  account)  since  the  contact 
angle  in  the  experiment  is  relatively  large  ($e  ~  70°). 

4.  Impact  of  wettability  on  evaporation  pattern 

Suppose  that  the  GDL  is  completely  saturated  by  water  after 
flooding.  One  obvious  option  to  restore  the  GDL  is  to  circulate  a 
dry  gas  within  the  bipolar  channel  so  as  to  dry  out  the  GDL.  We 
wish  to  explore  the  impact  of  wettability  on  drying  pattern  and 
drying  time.  To  this  end,  we  again  consider  2D  model  porous 
media  and  perform  drying  experiments  with  two  physical  mod¬ 
els  of  same  microstructure  but  of  different  wettability.  One  is 
machined  in  a  PTFE  plate  and  is  hydrophobic.  The  second  is 
machined  in  a  Plexiglas  plate  and  is  hydrophilic  (0Water  ~  80°). 
A  1-cm  thick  flat  Plexiglas  plate  is  fixed  on  top  of  the  machined 
plate.  This  transparent  plate  allows  for  the  direct  visualisation 
of  evaporation  pattern.  For  the  hydrophobic  system,  a  very  thin 
silicon  oil  film  carpets  the  top  Plexiglas  plate  so  as  to  make  it 
hydrophobic  (with  a  static  contact  angle  close  to  110°). 

In  continuity  with  previous  studies  on  drying  of  network 
micromodels,  we  do  not  adopt  a  disk  array  structure  as  before, 
but  rather  a  square  geometry  like  in  Ref.  [26].  This  does 
not  affect,  however,  the  findings  reported  in  this  section.  The 
microstructure  of  the  machined  models  is  shown  in  Fig.  10.  Note 
that  the  1-mm  depth  of  machined  ducts  is  uniform,  lower  than 
the  capillary  length  and  greater  than  the  duct  widths,  so  that  the 
invasion  is  essentially  controlled  by  the  duct  width  and  not  by 
its  depth.  The  channel  width  is  distributed  randomly  according 
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Fig.  10.  Model  porous  medium  used  for  the  evaporation  experiment. 
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Fig.  11.  Evaporation  patterns  (gas  phase  in  grey  or  white  (c),  water  in  black  or  dark  grey,  the  little  rectangles  of  random  size  form  the  solid  matrix),  (a)  Simulation 
of  evaporation  in  the  hydrophilic  network,  (b)  experimental  pattern  in  hydrophilic  network,  (c)  experimental  pattern  in  hydrophobic  network,  and  (d)  simulation  of 
evaporation  in  the  hydrophobic  network. 


to  a  discrete  uniform  distribution.  Except  for  the  uncertainties 
associated  with  the  machining  process,  five  classes  of  channel 
width  are  considered:  0.3,  0.4,  0.5,  0.6  and  0.7  mm.  A  5  mm 
in  width,  3  mm  in  depth  and  45  mm  long  channel  is  machined 
along  the  top  edge  of  the  network  as  shown  in  Fig.  10. 

The  network  is  initially  filled  with  distilled  water.  Evapora¬ 
tion  is  driven  by  the  flow  of  air  at  controlled  relative  humidity  in 
the  channel.  The  relative  humidity  (RH)  is  controlled  by  circulat¬ 
ing  the  air  in  a  chamber  containing  a  LiCl  solution  (RH  =  12%). 
The  system,  including  the  network  and  the  channel,  is  placed  in 
a  chamber  of  controlled  temperature  (40  °C).  The  experiments 
are  recorded  with  a  CCD  camera  set  above  the  chamber. 


The  obtained  evaporation  patterns  are  shown  in  Fig.  11. 
As  can  be  seen  from  Fig.  lib,  the  evaporation  pattern  in  the 
hydrophilic  model  resembles  a  capillary  fingering  pattern.  This 
result  was  expected.  From  previous  studies,  e.g.  Ref.  [27]  and 
references  therein,  it  is  known  that  slow  evaporation  and  quasi¬ 
static  drainage  lead  to  the  same  invasion  pattern  except  for  the 
isolated  clusters  that  are  trapped  for  good  in  drainage  but  not  in 
the  evaporation  situation  since  they  are  subject  to  evaporation. 
Fig.  1  la  shows  the  pattern  obtained  with  the  pore-network  evap¬ 
oration  model  proposed  in  Ref.  [28].  As  described  in  Ref.  [29], 
this  model  combines  the  invasion  percolation  algorithm  with 
the  computation  of  diffusion  transport  of  the  water  vapour  in  the 
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Fig.  12.  Sketch  of  liquid  invasion  and  O2  transport  in  a  hydrophilic  GDL  and  a  hydrophobic  one. 
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gas  phase.  As  can  be  seen  from  Fig.  11a  and  b,  the  agreement  is 
good  between  the  experimental  and  the  simulated  patterns.  As 
previously,  the  slight  differences  between  the  two  patterns  are 
attributed  to  the  effect  of  contact  angle,  which  is  on  the  order 
of  80°  in  the  experiment.  Again,  from  Section  3,  we  know  that, 
contrary  to  what  it  is  assumed  in  the  simulation,  the  invasion 
is  not  controlled  only  by  the  burst  mechanism  for  such  a  high 
value  of  contact  angle  (see  Figs.  8  and  9). 

Evaporation  in  the  hydrophobic  model,  see  Fig.  11c,  leads  to 
a  compact  pattern  with  an  almost  flat  invasion  front  reminiscent 
of  the  pattern  obtained  numerically  for  0e  >  100°  in  Section  3, 
see  Fig.  8.  For  the  hydrophobic  network,  air  is  the  wetting  fluid 
and  liquid  water  is  the  non- wetting  one.  Hence  slow  evapora¬ 
tion  in  this  case  is  analogous  to  the  process  of  displacement  of  a 
non- wetting  fluid  by  a  wetting  one.  This  process  is  usually  called 
imbibition,  e.g.  Ref.  [1 1].  It  is  therefore  not  surprising  that  evap¬ 
oration  leads  in  this  case  to  an  imbibition  pattern.  As  illustrated 
in  Fig.  lid,  it  is  possible  to  develop  pore-network  simulation 
of  evaporation  in  the  hydrophobic  network  by  combining  local 
imbibition  rules  and  the  computation  of  the  diffusion  transport 
in  the  gas  phase,  see  Ref.  [29]. 

The  great  difference  in  the  patterns  between  the  hydrophobic 
and  the  hydrophilic  networks  are  associated  with  differences  in 
drying  time  (for  the  same  condition  in  the  channel).  This  was 
explored  in  detail  in  Ref.  [29],  where  it  is  shown  that  drying  is 
faster  in  the  hydrophilic  network.  In  square  networks,  drying  can 
be  about  up  to  46%  faster  in  small  10x10  hydrophilic  networks 
compared  to  a  hydrophobic  one  of  the  same  size.  The  discrep¬ 
ancy  in  drying  time  decreases,  however,  with  the  network  size, 
with  drying  about  25%  faster  in  a  hydrophilic  network  for  large 
networks. 

5.  Discussion  and  conclusion 

As  regards  two-phase  flows  and  evaporation  in  GDLs,  the 
main  outcomes  of  the  results  presented  in  this  paper  can  be 
expressed  as  follows: 

•  Rendering  the  GDL  hydrophobic  is  a  good  option.  As 
sketched  in  Fig.  12,  large  regions  of  the  pore  space  remains 
accessible  to  the  gas  phase  up  to  the  CL/GDL  interface  in 
the  capillary  fingering  regime  that  characterizes  water  inva¬ 
sion  in  a  hydrophobic  GDL.  In  contrast,  the  compact  invasion 
in  a  hydrophilic  GDL  is  very  detrimental  to  the  gas  access 
to  the  GDL/CA  interface.  As  a  result,  the  liquid  saturation 
S  at  breakthrough  (when  the  liquid  reaches  the  channel)  in 
the  hydrophobic  GDL  is  much  lower  than  for  a  hydrophilic 
one  (where  Set  ~  1).  Thus  our  results  provide  a  pore-scale 
physical  explanation  to  the  fact  that  it  is  advantageous  to  use 
hydrophobic  GDLs. 

•  Drying  is  faster  with  a  hydrophilic  GDL.  However,  the  differ¬ 
ence  in  drying  time  with  a  hydrophobic  GDL  is  not  sufficient 
to  counterbalance  the  advantages  of  a  hydrophobic  GDL  in 
terms  of  liquid  distribution  during  the  GDL  invasion  by  water. 

•  The  capillary  number  characterizing  the  liquid  flow  in  a  GDL 
is  very  small.  Combined  with  the  fact  that  a  GDL  is  a  very  thin 
system,  the  flow  can  then  be  considered  as  dominated  by  cap¬ 


illary  effects.  To  simulate  the  liquid  flow  in  this  quasi-static 
regime  in  a  hydrophobic  GDL  at  the  scale  of  the  microstruc¬ 
ture,  one  can  rely  on  the  IP  algorithm  as  a  first  approximation. 
However,  as  shown  in  Section  3,  the  local  invasion  mecha¬ 
nisms  are  not  dominated  by  the  burst  mechanism  for  “large” 
contact  angles  0e  approaching  90° .  This  leads  to  changes  in 
the  invasion  pattern  compared  to  the  IP  pattern.  We  recall  that 
0e  ~  20°  for  the  air/water/PTFE  system  and  even  closer  to  90° 
within  a  GDL  material  according  to  [9] .  This  should  be  kept 
in  mind  when  comparing  experimental  data  and  simulations 
based  on  the  IP  algorithm  (or  the  algorithm  used  for  example 
in  Ref.  [15],  which  is  equivalent  to  IP). 

•  The  fact  that  the  invasion  regime  (compact  vs.  capillary  fin¬ 
gering)  is  completely  different  depending  on  the  wettability 
of  the  GDL  suggests  that  it  is  not  correct  to  use  basically  the 
same  Leverett  function  and  the  same  expressions  of  relative 
permeabilities  whatever  the  wettability  of  the  GDL,  e.g.  Ref. 
[16]. 

In  our  study,  only  2D  systems  have  been  considered.  It  is 
naturally  desirable  to  consider  3D  systems,  e.g.  Refs.  [13,15]. 
This  makes  significantly  more  difficult  the  numerical  simula¬ 
tions  and  the  visualisations.  The  2D  results  presented  here  may 
serve  as  guidelines  in  the  analysis  of  3D  results.  In  particular, 
all  the  results  listed  above  are  expected  to  be  qualitatively  valid 
in  3D. 

Since  a  MPL  is  often  attached  to  the  GDL,  it  is  interest¬ 
ing  to  study  two-phase  flows  in  MPL.  Using  multi- scale  pore 
networks  would  be  interesting  in  this  perspective.  In  the  same 
spirit,  pore-network  simulations  would  be  certainly  useful  to 
understand  better  the  complex  interactions  between  the  electro¬ 
chemical  and  transport  phenomena,  including  two-phase  flows 
and  wettability  effects,  in  catalyst  layers  as  well  as  the  interac¬ 
tions  between  the  various  layers  forming  a  PEMFC. 
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